\defmodule{NegativeMultinomialDist}

Implements the class \class{DiscreteDistributionIntMulti} for the
{\em negative multinomial} distribution with parameters $n > 0$ and
($p_1, \ldots, p_d$)  such that all $0<p_i<1$ and  $\sum_{i=1}^{d} p_i < 1$.
The probability mass function is \cite{tJOH69a}
\begin{htmlonly}
\eq
   P[X = (x_1, \ldots ,x_d)] = \left({\Gamma\left(n + \sum_{i=1}^{d} x_i\right)}/
        {\Gamma(n)}\right)
      p_0^{n}\prod_{i=1}^{d} p_i^{x_i}/x_i!
\endeq
\end{htmlonly}
\begin{latexonly}
\eq
   P[X = (x_1, \ldots,x_d)] = \frac{\Gamma\left(n +
                 \sum_{i=1}^{d} x_i\right)  p_0^{n}}
                 {\Gamma(n)}
    \prod_{i=1}^{d} \frac{p_i^{x_i}}{x_i!}
\eqlabel{eq:fNegativeMultinomial}
\endeq
\end{latexonly}
where $p_0 = 1 - \sum_{i=1}^{d} p_i$.

\bigskip\hrule

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{code}
\begin{hide}
/*
 * Class:        NegativeMultinomialDist
 * Description:  negative multinomial distribution
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.probdistmulti;
\begin{hide}
import umontreal.iro.lecuyer.util.Num;
import umontreal.iro.lecuyer.util.RootFinder;
import umontreal.iro.lecuyer.functions.MathFunction;
\end{hide}

public class NegativeMultinomialDist extends DiscreteDistributionIntMulti \begin{hide} {
   protected double n;
   protected double p[];

   private static class Function implements MathFunction {
      protected double Fl[];
      protected int ups[];
      protected int k;
      protected int M;
      protected int sumUps;

      public Function (int k, int m, int ups[], double Fl[]) {
         this.k = k;
         this.M = m;

         this.Fl = new double[Fl.length];
         System.arraycopy (Fl, 0, this.Fl, 0, Fl.length);
         this.ups = new int[ups.length];
         System.arraycopy (ups, 0, this.ups, 0, ups.length);

         sumUps = 0;
         for (int i = 0; i < ups.length; i++)
            sumUps += ups[i];
      }

      public double evaluate (double gamma) {
         double sum = 0.0;
         for (int l = 0; l < M; l++)
            sum += (Fl[l] / (gamma + (double) l));
         return (sum - Math.log1p (sumUps / (k * gamma)));
      }
   }


   private static class FuncInv extends Function implements MathFunction {

      public FuncInv (int k, int m, int ups[], double Fl[]) {
         super (k, m, ups, Fl);
      }

      public double evaluate (double nu) {
         double sum = 0.0;
         for (int l = 0; l < M; l++)
            sum += Fl[l] / (1.0 + nu * l);
         return (sum *nu - Math.log1p (sumUps * nu / k));
      }
   }
\end{hide}
\end{code}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Constructors}

\begin{code}

   public NegativeMultinomialDist (double n, double p[]) \begin{hide} {
      setParams (n, p);
   }\end{hide}
\end{code}
\begin{tabb}
   Creates a \texttt{NegativeMultinomialDist} object with parameters $n$
    and ($p_1$, \ldots, $p_d$) such that $\sum_{i=1}^{d} p_i < 1$,
   as described above. We have $p_i = $ \texttt{p[i-1]}.
\end{tabb}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Methods}

\begin{code}\begin{hide}

   public double prob (int x[]) {
      return prob_ (n, p, x);
   }
/*
   public double cdf (int x[]) {
      throw new UnsupportedOperationException ("cdf not implemented");
   }
*/
   public double[] getMean() {
      return getMean_ (n, p);
   }

   public double[][] getCovariance() {
      return getCovariance_ (n, p);
   }

   public double[][] getCorrelation() {
      return getCorrelation_ (n, p);
   }

   private static void verifParam (double n, double p[]) {
      double sumPi = 0.0;

      if (n <= 0.0)
         throw new IllegalArgumentException("n <= 0");

      for (int i = 0; i < p.length;i++) {
         if ((p[i] < 0) || (p[i] >= 1))
            throw new IllegalArgumentException("p is not a probability vector");

         sumPi += p[i];
      }
      if (sumPi >= 1.0)
         throw new IllegalArgumentException("p is not a probability vector");
   }

   private static double prob_ (double n, double p[], int x[]) {
      double p0 = 0.0;
      double sumPi= 0.0;
      double sumXi= 0.0;
      double sumLnXiFact = 0.0;
      double sumXiLnPi = 0.0;

      if (x.length != p.length)
         throw new IllegalArgumentException ("x and p must have the same size");

      for (int i = 0; i < p.length;i++)
      {
         sumPi += p[i];
         sumXi += x[i];
         sumLnXiFact += Num.lnFactorial (x[i]);
         sumXiLnPi += x[i] * Math.log (p[i]);
      }
      p0 = 1.0 - sumPi;

      return Math.exp (Num.lnGamma (n + sumXi) - (Num.lnGamma (n) +
           sumLnXiFact) + n * Math.log (p0) + sumXiLnPi);
   }\end{hide}

   public static double prob (double n, double p[], int x[])\begin{hide} {
      verifParam (n, p);
      return prob_ (n, p, x);
   }\end{hide}
\end{code}
\begin{tabb}
   Computes the probability mass function (\ref{eq:fNegativeMultinomial})
   of the negative multinomial distribution with parameters
   $n$ and ($p_1$, \ldots, $p_d$), evaluated at $\mathbf{x}$.
\end{tabb}
\begin{code}\begin{hide}

   private static double cdf_ (double n, double p[], int x[]) {
      throw new UnsupportedOperationException ("cdf not implemented");
   }\end{hide}

   public static double cdf (double n, double p[], int x[])\begin{hide} {
      verifParam (n, p);

      return cdf_ (n, p, x);
   }\end{hide}
\end{code}
\begin{tabb}
   Computes the cumulative probability function $F$ of the
   negative multinomial distribution with parameters $n$
   and ($p_1$, \ldots, $p_k$), evaluated at $\mathbf{x}$.
\end{tabb}
\begin{code}\begin{hide}

   private static double[] getMean_ (double n, double p[]) {
      double p0 = 0.0;
      double sumPi= 0.0;
      double mean[] = new double[p.length];

      for (int i = 0; i < p.length;i++)
         sumPi += p[i];
      p0 = 1.0 - sumPi;

      for (int i = 0; i < p.length; i++)
         mean[i] = n * p[i] / p0;

      return mean;
   }\end{hide}

   public static double[] getMean (double n, double p[])\begin{hide} {
      verifParam (n, p);

      return getMean_ (n, p);
   }\end{hide}
\end{code}
\begin{tabb}
   Computes the mean $E[X] = n p_i / p_0$ of the negative multinomial distribution
   with parameters $n$ and ($p_1$, \ldots, $p_d$).
\end{tabb}
\begin{code}\begin{hide}

   private static double[][] getCovariance_ (double n, double p[]) {
      double p0 = 0.0;
      double sumPi= 0.0;
      double cov[][] = new double[p.length][p.length];

      for (int i = 0; i < p.length;i++)
         sumPi += p[i];
      p0 = 1.0 - sumPi;

      for (int i = 0; i < p.length; i++)
      {
         for (int j = 0; j < p.length; j++)
            cov[i][j] = n * p[i] * p[j] / (p0 * p0);

         cov[i][i] = n * p[i] * (p[i] + p0) / (p0 * p0);
      }

      return cov;
   }\end{hide}

   public static double[][] getCovariance (double n, double p[])\begin{hide} {
      verifParam (n, p);

      return getCovariance_ (n, p);
   }\end{hide}
\end{code}
\begin{tabb}
   Computes the covariance matrix of the negative multinomial distribution
   with parameters $n$ and ($p_1$, \ldots, $p_d$).
\end{tabb}
\begin{code}\begin{hide}

   private static double[][] getCorrelation_ (double n, double[] p) {
      double corr[][] = new double[p.length][p.length];
      double sumPi= 0.0;
      double p0;

      for (int i = 0; i < p.length;i++)
         sumPi += p[i];
      p0 = 1.0 - sumPi;

      for (int i = 0; i < p.length; i++) {
         for (int j = 0; j < p.length; j++)
            corr[i][j] = Math.sqrt(p[i] * p[j] /((p0 + p[i]) * (p0 + p[j])));
         corr[i][i] = 1.0;
      }
      return corr;
   }\end{hide}

   public static double[][] getCorrelation (double n, double[] p)\begin{hide} {
      verifParam (n, p);
      return getCorrelation_ (n, p);
   }\end{hide}
\end{code}
\begin{tabb}
   Computes the correlation matrix of the negative multinomial distribution
   with parameters $n$ and ($p_1$, \ldots, $p_d$).
\end{tabb}
\begin{code}

   public static double[] getMLE (int x[][], int m, int d)\begin{hide} {
      int ups[] = new int[m];
      double mean[] = new double[d];

      int i, j, l;
      int M;
      int prop;

      // Initialization
      for (i = 0; i < d; i++)
         mean[i] = 0;

      // Ups_j = Sum_k x_ji
      // mean_i = Sum_m x_ji / m
      for (j = 0; j < m; j++) {
         ups[j] = 0;
         for (i = 0; i < d; i++) {
            ups[j] += x[j][i];
            mean[i] += x[j][i];
         }
      }
      for (i = 0; i < d; i++)
         mean[i] /= m;

/*
      double var = 0.0;
      if (d > 1) {
         // Calcule la covariance 0,1
         for (j = 0; j < m; j++)
            var += (x[j][0] - mean[0])*(x[j][1] - mean[1]);
         var /= m;
      } else {
         // Calcule la variance 0
         for (j = 0; j < m; j++)
            var += (x[j][0] - mean[0])*(x[j][0] - mean[0]);
         var /= m;
      }
*/

      // M = Max(Ups_j)
      M = ups[0];
      for (j = 1; j < m; j++)
         if (ups[j] > M)
            M = ups[j];

      if (M >= Integer.MAX_VALUE)
         throw new IllegalArgumentException("n/p_i too large");

      double Fl[] = new double[M];
      for (l = 0; l < M; l++) {
         prop = 0;
         for (j = 0; j < m; j++)
            if (ups[j] > l)
               prop++;

         Fl[l] = (double) prop / (double) m;
      }

/*
      // Estime la valeur initiale de n pour brentDekker: pourrait
      // accélérer brentDekker (gam0/1000, gam0*1000, f, 1e-5).
      // Reste à bien tester.
      if (d > 1) {
         double gam0 = mean[0] * mean[1] / var;
         System.out.println ("gam0 = " + gam0);
      } else {
         double t = var/mean[0] - 1.0;
         double gam0 = mean[0] / t;
         System.out.println ("gam0 = " + gam0);
      }
*/
      double parameters[] = new double[d + 1];
      Function f = new Function (m, (int)M, ups, Fl);
      parameters[0] = RootFinder.brentDekker (1e-9, 1e9, f, 1e-5);

      double lambda[] = new double[d];
      double sumLambda = 0.0;
      for (i = 0; i < d; i++) {
         lambda[i] = mean[i] / parameters[0];
         sumLambda += lambda[i];
      }

      for (i = 0; i < d; i++) {
         parameters[i + 1] = lambda[i] / (1.0 + sumLambda);
         if (parameters[i + 1] > 1.0)
            throw new IllegalArgumentException("p_i > 1");
      }

      return parameters;
   }\end{hide}
\end{code}
\begin{tabb}
   Estimates  and returns the parameters [$\hat{n}$, $\hat{p}_1$, \ldots,
  $\hat{p}_d$]
   of the negative multinomial distribution using the maximum likelihood method.
   It uses the $m$ observations of $d$ components in table
   \texttt{x[$i$][$j$]}, $i = 0, 1, \ldots, m-1$ and $j = 0, 1, \ldots, d-1$.
   \begin{detailed}
   The equations of the maximum likelihood are defined in \cite{tJOH69a}:
   \begin{eqnarray*}
      \sum_{s=1}^{M} \frac{F_s}{(\hat{n} + s - 1)} & = & \ln \left(1 +
        \frac{1}{\hat{n} m} \sum_{j=1}^{m} \Upsilon_j \right)\\[8pt]
      p_i & = & \frac{\lambda_i}{1 + \sum_{j=1}^{d} \lambda_j}
                 \qquad \mbox{for } i=1, \ldots,d
   \end{eqnarray*}
   where
   \begin{eqnarray*}
      \lambda_i & = & \frac{\sum_{j=1}^{m} X_{i,j}}{\hat{n} m}
             \qquad \mbox{for } i=1, \ldots,d\\
      \Upsilon_j & = & \sum_{i=1}^{d} X_{i,j} \qquad \mbox{for } j=1, \ldots,m\\
      F_s & = & \frac{1}{m} \sum_{j=1}^{m} \mbox{\textbf{1}}\{\Upsilon_j \ge s\}
        \qquad \mbox{for } s=1, \ldots,M\\
      M & = & \max_j \{\Upsilon_j\}
   \end{eqnarray*}
   \end{detailed}
\end{tabb}
\begin{htmlonly}
   \param{x}{the list of observations used to evaluate parameters}
   \param{m}{the number of observations used to evaluate parameters}
   \param{d}{the dimension of each vector}
   \return{returns the parameters [$\hat{n}$, $\hat{p}_1$, \ldots, $\hat{p}_d$]}
\end{htmlonly}
\begin{code}

   public static double getMLEninv (int x[][], int m, int d)\begin{hide} {
      int ups[] = new int[m];
      double mean[] = new double[d];
      int i, j, l;
      int M;
      int prop;

      // Initialization
      for (i = 0; i < d; i++)
         mean[i] = 0;

      // Ups_j = Sum_k x_ji
      // mean_i = Sum_m x_ji / m
      for (j = 0; j < m; j++) {
         ups[j] = 0;
         for (i = 0; i < d; i++) {
            ups[j] += x[j][i];
            mean[i] += x[j][i];
         }
      }
      for (i = 0; i < d; i++)
         mean[i] /= m;

      // M = Max(Ups_j)
      M = ups[0];
      for (j = 1; j < m; j++)
         if (ups[j] > M)
            M = ups[j];

      if (M >= Integer.MAX_VALUE)
         throw new IllegalArgumentException("n/p_i too large");

      double Fl[] = new double[M];
      for (l = 0; l < M; l++) {
         prop = 0;
         for (j = 0; j < m; j++)
            if (ups[j] > l)
               prop++;

         Fl[l] = (double) prop / (double) m;
      }

      FuncInv f = new FuncInv (m, M, ups, Fl);
      double nu = RootFinder.brentDekker (1.0e-8, 1.0e8, f, 1e-5);
      return nu;
   }\end{hide}
\end{code}
\begin{tabb}
   Estimates  and returns the parameter $\nu = 1/\hat{n}$
   of the negative multinomial distribution using the maximum likelihood method.
   It uses the $m$ observations of $d$ components in table
   \texttt{x[$i$][$j$]}, $i = 0, 1, \ldots, m-1$ and $j = 0, 1, \ldots, d-1$.
   \begin{detailed}
   The equation of maximum likelihood is defined as:
   \[
      \sum_{s=1}^{M} \frac{\nu F_s}{(1 + \nu(s - 1))}  =  \ln \left(1 +
        \frac{\nu}{m} \sum_{j=1}^{m} \Upsilon_j \right)\\[8pt]
   \]
   where the symbols are defined as in \method{getMLE}{int[][], int, int}{}.
   \end{detailed}
\end{tabb}
\begin{htmlonly}
   \param{x}{the list of observations used to evaluate parameters}
   \param{m}{the number of observations used to evaluate parameters}
   \param{d}{the dimension of each vector}
   \return{returns the parameter $1/\hat{n}$}
\end{htmlonly}
\begin{code}

   public double getGamma()\begin{hide} {
      return n;
   }\end{hide}
\end{code}
\begin{tabb}
   Returns the parameter $n$ of this object.
\end{tabb}
\begin{code}

   public double[] getP()\begin{hide} {
      return p;
   }\end{hide}
\end{code}
\begin{tabb}
   Returns the parameters ($p_1$, \ldots, $p_d$) of this object.
\end{tabb}
\begin{code}

   public void setParams (double n, double p[])\begin{hide} {
      if (n <= 0.0)
         throw new IllegalArgumentException("n <= 0");

      this.n = n;
      this.dimension = p.length;
      this.p = new double[dimension];

      double sumPi = 0.0;
      for (int i = 0; i < dimension; i++) {
         if ((p[i] < 0) || (p[i] >= 1))
            throw new IllegalArgumentException("p is not a probability vector");

         sumPi += p[i];
         this.p[i] = p[i];
      }

      if (sumPi >= 1.0)
         throw new IllegalArgumentException("p is not a probability vector");
   }\end{hide}
\end{code}
\begin{tabb}
   Sets the parameters $n$ and ($p_1$, \ldots, $p_d$) of this object.
\end{tabb}
\begin{code}\begin{hide}
}\end{hide}
\end{code}
